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The semimetal to antiferromagnet quantum phase transition of the Hubbard model on the honey¬ 
comb lattice has come to the forefront in the context of the proposal that a semimetal to spin liquid 
transition can occur before the transition to the antiferromagnetic phase. To study the semimetal 
to antiferromagnet transition, we generalize the two-particle self-consistent (TPSC) approach to the 
honeycomb lattice (a structure that can be realized in graphene for example). We show that the 
critical interaction strength where the transition occurs is Uc/t — 3.79 ±0.01 quite close to the value 
Celt — 3.869 ±0.013 reported using large-scale quantum Monte Carlo simulations. This reinforces 
the conclusion that the semimetal to spin liquid transition is pre-empted by the transition to the an¬ 
tiferromagnet. Since TPSC satisfies the Mermin-Wagner theorem, we find temperature-dependent 
results for the antiferromagnetic and ferromagnetic correlation lengths as well as the dependence 
of double occupancy and of the renormalized spin and charge interactions on the bare interaction 
strength. We also estimate the value of the crossover temperature to the renormalized classical 
regime as a function of interaction strength. 

I. INTRODUCTION 

Consider a half-filled band of electrons that interact 
through a short-range potential 7/ on a lattice with band¬ 
width W. As one increases interactions, the ground state 
can undergo a transition from a Fermi liquid to a mag¬ 
netically ordered state for U less than W, but if there is 
enough frustration, quantum fluctuations may prohibit 
long-range order. In that case, upon increasing U fur¬ 
ther there may be a Fermi liquid to insulator transition 
(a Mott transition) where the insulator, a spin liquid, 
does not exhibit long-range order. 

Spin liquids have been extensively searched for since 
Anderson’s proposal in the context of high-temperature 
superconductors The pyrochlore spin ices with frac¬ 
tionalized excitations are the best candidates to date for 
spin liquid ground states in three dimensions Since 
quantum fluctuations are large in low dimension, two- 
dimensional lattices are especially good candidates for 
spin liquid ground states. There is experimental evidence 
for such a state of matter in layered organic materials of 
the BEDT family that form a highly frustrated triangular 
lattice The first theoretical proposal for a spin liquid 
state in the 1970’s was in fact for the triangular lattice 

Theoretically, evidence for a spin liquid state has also 
been found on the kagome lattice 

The honeycomb lattice stands as a particularly in¬ 
teresting candidate for a spin liquid ground state be¬ 
cause it has the smallest possible coordination for a two- 


dimensional lattice, leading to large quantum fluctua¬ 
tions. In addition, the Hubbard model on the honeycomb 
lattice may be relevant for a number of real systems, in¬ 
cluding graphene, carbon nanotubes, MgB 2 , etc., as men¬ 
tioned in Ref. 6. 

Much recent work has focused on this model ever since 
very large scale quantum Monte Carlo simulations made 
the exciting prediction of a spin liquid over a small range 
of values of U, beyond which antiferromagnetism sets 
in.^ This claim has been confirmed by further numerical 
work^“^^ but was later disputed by Sorella et using 
even larger lattices. 

Since methods based on dynamical mean-field theory 
(DMFT) and its extensions^^~^^-so-called quantum clus¬ 
ter approaches-are particularly suited to find Mott tran¬ 
sitions, they have been used to look for a spin liquid phase 
between the semimetal and the antiferromagnet. After 
early single-site DMFT studies quantum cluster 

calculations confirmed the existence of the intermediate 
spin iiquid phase^^"^^ or of a Mott transition^^. How¬ 
ever, Hassan et using the ciuster dynamical im¬ 

purity approximation (GDIA) found out that the Mott 
transition necessary for a spin liquid ground state is in 
fact pre-empted by antiferromagnetic long-range order. 
Careful analysis^^’^^ of the influence of the cluster shape 
and of the various implementations of cluster extensions 
of dynamical mean-field theory^^’^^ suggest that it is im¬ 
portant to apply other quantitative methods to find the 
precise values of the critical values of UjW for the phase 
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transitions^^. Other approaches that have been applied 
to this problem are briefly summarized in Refs. 25 and 
29. 

Since the Mott transition towards a spin liquid occurs 
in the absence of long-range order, one can expect that 
quantum cluster methods give a good upper bound for 
the occurrence of this transition^^. However, to locate 
the precise value of Uc/W where antiferromagnetism sets 
in, it is crucial that the method correctly include long- 
wavelength quantum fluctuations in the thermodynamic 
limit. Given that the critical value of UdW for antifer¬ 
romagnetism is of order 2/3 a semianalytical, nonper- 
turbative technique, valid from weak to intermediate cou¬ 
pling, the two-particle self consistent approach, (TPSC) 
is especially suited for this problem^^’^^. This is the ap¬ 
proach we use in this paper. Unlike RPA or Hartree- 
Fock theory, this method satisfies not only conservation 
laws, but also the Pauli principle, the Mermin-Wagner 
theorem, and important sum rules for spin and charge 
fluctuations. TPSC allows us to locate the crossover 
to the renormalized classical regime where the correla¬ 
tion length for antiferromagnetic fluctuations exceeds the 
thermal de Broglie wavelength. The extrapolation of that 
crossover line to zero temperature is one of the methods 
that can be used to find the value of UdW where an¬ 
tiferromagnetism sets in. While TPSC has so far been 
used only in a single-band context, here we generalize it 
to the two-band case to find UdW. 

Previous estimates of Uc for the antiferromagnetic 
transition of the Hubbard model with a nearest-neighbor 
hopping t on the honeycomb lattice at half-filling and 
T = 0 are in the 3.5t to 5t, much 

larger than the Hartree-Fock RPA mean-field result 
2.23t. A number of numerical lattice-field theory solu¬ 
tions of the continuum problem (see Ref. 24 of Liebsch 
and Wu^^) also suggest an antiferromagnetic phase (more 
precisely chiral symmetry breaking) at strong coupling. 
The most accurate estimate for Uc should be the recent 
large scale quantum Monte Carlo calculation of Sorella 
et Uc/t = 3.869 ± 0.013, close to estimates from 

high temperature series expansion^, Uc ~ 4t and from a 
projection Monte Carlo method with an optimized initial 
state by Furukawa^^ Uc ^ 3.6t. Another accurate recent 
result, Uc = 3.78t, is provided by the pinning field ap¬ 
proach to quantum Monte Carlo of Assaad and Herbut 
Other quantum Monte Carlo calculations generally 
find higher values Uc. This includes the early ones by 
Sorella and Tosatti^^ that yielded Uc = 4.5t, those of 
Paiva et al.^ that found Uc ~ 5t and those of Meng et 
al.^ with Uc > 4.3t. The most accurate weak-coupling 
method that can be compared with TPSC, namely, the 
functional renormalization group^^’^^, gives Uc ~ 3.8t, 
close to the best estimates mentioned above. 

The paper is organized as follows. In Sec. H, we in¬ 
troduce the model and the notation for the Green func¬ 
tion formalism. We generalize the TPSC approach to 
graphene in Sec. HI, obtaining the spin and charge fluctu¬ 
ations with a functional derivative approach. The scaling 


for the susceptibility is obtained in Sec. IV. The numer¬ 
ical procedure is explained in Sec V and the numerical 
results are presented in Sec. VI. Three appendices con¬ 
tain analytical results that can be obtained for the spin 
susceptibility. 


II. MODEL AND GREEN FUNCTION 


The Hamiltonian is given by 

n = Ho + U'^ ni^jiii ( 1 ) 

i 

Ho = -t + ( 2 ) 

<ij>a 


where Hq is the noninteracting hopping Hamiltonian. 
Creation operators for a particle on sublattice A and B 
are represented by d and 6^, respectively, a is the spin 
of the particle and < ij > represents nearest-neighbor 
sites on the honeycomb lattice. Here, t is the hopping 
parameter and U is the strength of the on-site Coulomb 
interaction. 

In Fourier space, Hq takes the form 


where 


Hq = 



J 


/(k) = 1 + 


( 3 ) 

( 4 ) 


with ai = and a 2 = the basis vec- 

tors of length unity for the underlying triangular Bravais 
lattice. We take the nearest-neighbor hopping t equal to 
unity. Similarly, the Planck’s constant h and Boltzmann 
constant ks are set to unity. 

The Green function for the Hamiltonian in Eq. (I) is 
a 4 X 4 matrix since there are two sublattices and two 
spin indices. The Green function matrix G is diagonal in 
spin-space because of the spin rotational invariance of the 
Hamiltonian (I). Introducing the notation I = (ri,ri), 
where I stands for the position on the triangular lattice 
rl and imaginary time ri, the matrix elements of G are 
defined by 


GS^'(l,2) = (5) 

where a = ad and /3 = a, 6 denote sublattice indices and 
cr, a' =t, i are the spin indices. 

The equation of motion for (1? 2) in Eq. (5) is 


5G-5{1,2) 


dn 


(5(ti '^2)^rir2^crcr'^a/3 

{Tr-da,{l)pl{2))5,,,. 


( 6 ) 
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The Heisenberg equation of motion in the grand canoni¬ 
cal ensemble yields 

-^a„{l) = [H- (7) 

OTi 

where fi is the chemical potential and M is the total- 
number operator. Defining 

= -t^ <^ri+A,r2 5{ti - , (8) 

A 

where P = a^b are the sublattice indices, A runs over 
the nearest neighbors, and is the Pauli matrix 

c=(; 5). (s» 

the equation of motion for the Green function takes the 
form 

(-^ + M)4.r3GS^'(3,2) - (l,3)G^y(3,2) 

= S(^Ti '^2)^rir2^crcr'^a/3 

- [/(T,4(l)a<x(l)a<x(l)^i(2))(5 aa' 5 (10) 

where a bar over an index like 3 implies summation over 
the corresponding lattice positions and an integral over 
imaginary time, while the Einstein summation conven¬ 
tion applies to repeated spin or sublattice indices. 

Using 

Go\l,2) = {-dr,+fi)l-h{l,2), (11) 

for the noninteracting Green function, a short-hand for 
the above equation of motion is 

Go 1(1, 3)G(3, 2) = 5(n - T 2 )Sr,r,l - u(l, 2) (12) 

where the four-point correlation matrix u is 

<^'(1,2) = -f/(T,4(l+)a,,(l)a,,(l)^t(2))<5^^,. (13) 

The correlation matrix can be rewritten in terms of the 
self-energy using Dyson’s equation 

G-i(l,2) = Go'(l,2)-S(l,2) (14) 

where from the spin-symmetry of the Hamiltonian, the 
self-energy is block-diagonal in spin subspace. This leads 
to 

S(1,3)G(3,2) =u(l,2). (15) 

This clearly shows how the self-energy is related to two- 
particle correlation functions and to the potential energy 
in the special case where the first and last indices are 
equal (with a small positive shift in imaginary time for 
proper time-order). This well-known relation is obtained 
without any approximations. This is the multi-band gen¬ 
eralization of an important consistency requirement be¬ 
tween the self-energy and the double occupancy in the 
Hubbard model^^. 


III. GENERALIZATION OF TPSC 


The two-particle self-consistent (TPSC) approach 
was developed to study the single-band Hubbard 
model^^’^^’^^~^^. It has been benchmarked through de¬ 
tailed comparisons with quantum Monte Carlo calcula¬ 
tions. This is a nonperturbative method that works best 
from weak to intermediate values of coupling UjW. The 
key features of this approach are that it satisfies conser¬ 
vation laws, the Pauli principle and the Mermin-Wagner 
theorem. 

Perturbative methods, which obey conservation laws 
(like FLEX^^), tend to violate the Pauli principle, while 
those that satisfy the Pauli principle (like the Parquet re¬ 
summations) usually violate conservation laws^^. Meth¬ 
ods like RPA give a finite-temperature transition to an 
antiferromagnetic state with the long-range order, a sce¬ 
nario prevented by the Mermin-Wagner theorem in two 
dimensions 

Although the spin and charge susceptibilities in TPSC 
are similar in form to those appearing in RPA, the two 
methods are fundamentally different. In TPSC, the irre¬ 
ducible spin and charge vertices are not equal. They are 
assumed to be momentum and frequency independent 
and are computed self-consistently at the two-particle 
level in such a way that local sum rules for spin and 
charge are satisfied. With TPSC, one can study the an¬ 
tiferromagnetic fluctuations in two-dimensional lattices 
without the unphysical finite-temperature phase transi¬ 
tion. The renormalized classical regime, where the fluc¬ 
tuations are large and the correlation length becomes 
greater than the thermal de Broglie wavelength, can be 
studied using this theory. This crossover to a renormal¬ 
ized classical regime at a finite temperature is a precursor 
of the zero-temperature instability to the long-range or¬ 
der. Note, however, that TPSC is not valid deep inside 
the renormalized classical regime. 

TPSC has been used to demonstrate, for example, that 
antiferromagnetic fluctuations can induce a pseudogap 
in two dimensions"^^’^^’^^ and that d-wave superconduc¬ 
tivity mediated by these fluctuations is possible^^. The 
method has been generalized to the attractive Hubbard 
model,and has been extended to the case where one 
includes a near-neighbor repulsion V. This is called ex¬ 
tended TPSC, or ETPSC^^-^^ 

Here we generalize the method to two bands, but iden¬ 
tical atoms within the unit cell. We do not consider the 
second step of the theory, which gives an improved for¬ 
mula for the self-energy^^. The final form of the theory 
is very natural but we give a detailed derivation below. 
The reader may skip to the results section without loss 
of continuity. The relevant equations are the spin and 
charge susceptibilities (26) and (27) and the local spin 
and charge sum rules (29) and (30) that have to be solved 
self-consistently with the ansatz (18) and (23). 
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A. TP SC ansatz for two bands 

The renormalized interactions for spin and charge can 
be obtained from functional derivatives of the self-energy 
y; 3 i, 4 i obtain 5 ] from its definition in terms of 

the four-point function u in Eq. (15), we assume that 
Hartree-Fock factorization as a product of two-point cor¬ 
relation functions is justified when the four points in u do 
not coincide But when all four points in u, Eq. (13), 
are identical, we impose the exact relation given by 

[S(l,3)G(3,l + )];^^' SafsSa.' =U{naail)naM))Sa0S,.', 

(16) 

obtained from Eq. (15) when 2^1+ and a = /3, i.e. the 
positions coincide and the times are such that r 2 = ti -he, 
where e is positive and infinitesimal. 

Using spin-rotational invariance, since (16) is diago¬ 
nal in spin indices, we focus on the diagonal elements 
and then, the above Hartree-Fock-like factorization of 
Eq. (15) can be written as 

3)G;^(3, 2) = 2). (17) 

For 2^1+ and a = P, we find that the exact result 
(16) is recovered if 


This expression for A involves double occupancy 
{riaa'^aa)^ which is obtained by the self-consistent cal¬ 
culations explained in the next subsection. 

Substituting for A in Eq.(17) and right multiplying by 
we obtain 


ES^(l, 2) = 1+)<5(1 - 2)5a0. (19) 

This is our first approximation for the self-energy. It is lo¬ 
cal and frequency independent. A better approximation 
can be obtained by including the effects of fluctuations 
but this is not needed here^^’^^’^^. As explained in the 
next section, the functional derivatives of the self-energy 
obtained above lead to the renormalized vertices for spin 
and charge. 


B. Spin and charge susceptibilities 

The spin and charge susceptibilities are calculated to 
reach an understanding of the competing spin and charge 
ordering transitions in the model. The value of A in 
Eq. (18) is obtained from these susceptibilities. Unlike 
RPA, where vertices are the bare U in both spin and 
charge susceptibilities, in TPSC^^’^^, spin and charge ver¬ 
tices differ. The renormalized irreducible vertex for spin 
is denoted by Us and for charge by Uc- They will clearly 
both depend on U. 


The spin and charge vertices in the longitudinal spin 
channel are computed from the local particle-hole irre¬ 
ducible vertices and These vertices are given 

by functional derivatives of the self-energy 


-pcrcr' 

^ a/3,7C 


(1,2; 3,4) 


<^G-'(3,4)- 


( 20 ) 


In matrix notation for the sublattice indices, the irre¬ 
ducible spin vertex is given by 


r.(l,2;3,4) 


(5Ef(l,2) (5Ef(l,2) 

(5Gf(3,4) “ (5Gf(3,4)’ 


( 21 ) 


where ajS corresponds to the row index and corre¬ 
sponds to the column index with a, 7 , (" = a, b. From 
our first approximation for the self-energy, Eq. (19), we 
can calculate these functional derivatives. We check that 
the functional derivatives of the terms cancel 

{ncxt){nai) 

by spin-rotational invariance, and we obtain 


r,(l, 2 ; 3, 4 ) = U, (5(1 - 3 ) 5 ( 1 + - 4 ) 5(1 - 2 ), ( 22 ) 


where the only non-zero elements of Us are diagonal in 
sublattice index and are given by 


iiOLOL.acx — A — TT \^<Aa'^aa) 

* “ {naa){naa)' 

The irreducible charge vertex is 


rc(l,2;3,4) 


5S'f'(l,2) 5St(l,2) 

5G1(3,4) 5G1(3,4)' 


(23) 


(24) 


From the functional derivative of the terms , 

{rioctJAcxi) ’ 

we obtain correlation functions of higher order. TPSC 
makes the assumption that the irreducible charge vertex, 
like the irreducible spin vertex, is constant and diagonal 
in sublattice index: 


r^l, 2 ; 3, 4)=UcS(l- 3)5(1+ - 4)5(1 - 2 ). (25) 

Introducing the short hand q = (q, iz^), which stands 
for the momentum space coordinate q and the bosonic 
Matsubara frequency z/, we find a straightforward gen¬ 
eralization of the particle-hole Bethe-Salpeter equation^^ 
to the case of a matrix susceptibility. The corresponding 
spin susceptibility the charge susceptibility 

are given by. 


X%q) = 

(l-^X°{qn'^ X^q), 

(26) 

X%q) = 

(l+^X%q)Gc) X°{q), 

(27) 


where x^ is the noninteracting susceptibility (Lindhard 
function) defined by 

Xa/3,7C(9) = (28) 

kaa' 
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The summation is over the momentum space k as well as 
over the fermionic Matsubara frequencies, and the lattice 
size is X 

The sum rules^^ needed for self-consistency are ob¬ 
tained by summing susceptibilities over all momenta and 
frequencies to recover local equal-time correlation func¬ 
tions. In the spin channel, we find 

Jp^Xaa,aa{<l) = («at) + Ki) “ ■ (29) 

q 

On the right-hand side, we have used the fact that the 
Pauli principle must be satisfied in the form = n^j. 
The corresponding sum rules for the charge susceptibility 
are 

Jp X'aa,aa{^) = {^1) + (^i) + 2«ni) - {Uaf. 

(30) 

We already have an expression, Eq. (23), for Us in 
terms of double occupancy. By substituting this in 
Eq. (26) for the spin susceptibility, we can evaluate the 
sum rules given by Eq. (29) and obtain the values of the 
double occupancies {ua^Ua]) and and hence 

in a self-consistent manner. By symmetry, here {na\na\) 
and are equal. We can determine the constant 

charge vertex Uc from the sum rules given by Eq. (30) 
once we know the values of double occupancies. Now 
that we have Ug and Uc^ the susceptibilities can be cal¬ 
culated from Eqs. (26) and (27). 


We can study the fluctuations in the system as a func¬ 
tion of temperature T and on-site interaction U. The 
correlation lengths corresponding to various channels in 
the spin and charge susceptibilities give us an estimate 
of the magnitudes of the fluctuations and hence let us 
determine which ordering transition is dominant in the 
system. The crossover to a renormalized classical regime 
at lower temperatures can be detected from the corre¬ 
sponding correlation length. 

IV. SCALING FORM FOR THE 

SUSCEPTIBILITIES 

The correlation length is useful to find the renormal¬ 
ized classical regime and the zero-temperature critical 
value of U. In the limit where the correlation lengths 
are large, a simple analytical form is useful. First, we 
introduce the notation 

Xaa,aa Xaa 5 Xaa,bb Xab Xbb,aa = xl- (31) 
Since the a and b sublattices are equivalent, we will use 

Xaa,aa — Xbb,bb — Xaa' (32) 

Quite generally, we also have the following equality 

{xLMY = (33) 

where is a bosonic Matsubara frequency. 

The spin susceptibilities can conveniently be rewritten 
in terms of susceptibilities that are either ferromagnetic 
or antiferromagnetic within a unit cell. First, rewrite the 
determinant entering the spin susceptibility (26) as 


det (l - ( 9 )^/ 5 ) = 


1 - “ \Jxlb(.q)xlai(l)) 1 - + \Jxlbi^)xlai(l)) 


(34) 


with an analogous result for the determinant entering the 
charge susceptibility. Clearly, the location of the poles is 
determined by the combinations of noninteracting sus¬ 
ceptibilities 

x)t = (^Xaa - \/XabXba^ , (35) 

Xafm = + \JXabXba^ ■ (36) 

These can be associated with the noninteracting ferro¬ 
magnetic and antiferromagnetic spin susceptibilities, re¬ 
spectively. Note that the usual definition of antiferro¬ 
magnetism, which we adopt here, corresponds to alter¬ 
nating spin directions on a and b sublattices but occurs 
at q = 0 as far as wave vectors are concerned. Explicit 
expressions for the noninteracting susceptibilities appear 


in Appendix A. Intraband terms contribute to the ferro¬ 
magnetic susceptibility while the antiferromagnetic sus¬ 
ceptibility involves interband transitions. 

Taking the analogous definition for the interacting case 
we find, after some algebra detailed in Appendix B, the 
following scalar equations 


x/m(q>*0 


Xa/m(q>*0 


x/m(q^^O 
1- ^x/„(q,*0 

1- ^Xa/™(q,*0' 


(37) 

(38) 


They resemble the expressions in the single-band case. 
Analogous definitions can be made for the charge sus¬ 
ceptibilities. 

The correlation length becomes large when the denom- 
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inator of the interacting susceptibilities is close to zero 
at vanishing Matsubara frequency. Taking the antiferro- 
magnetic susceptibility as an example, in that situation 
the numerator replaced by the max¬ 
imum value Xafmi^ = = 0) while in the denomina¬ 
tor iz/) must be expanded about the maximum 

at q = 0, = 0. 

Because of the Dirac cones, the noninteracting suscep¬ 
tibility in the denominator does not have a derivative at 
q = 0, iz^ = 0. The left derivative and the right derivative 
as we approach q = 0 are different. In Appendix C, we 
estimate the derivatives in the Dirac approximation. 

We can proceed numerically to confirm the orders 
of magnitude obtained in Appendix C. From the con¬ 
ical shape of the surface plot (see Fig. 1) of the anti¬ 
ferromagnetic susceptibility the dependence is on q = 

\JQx ^ Qy' The derivative is obtained by fitting the data 

for x^j^(q, 0) about q = 0, using a form given by 

Xa/m(q>0) =a + ^\jQi + ll+c{ql + ql)- The coefficient 

b = dXafml^^ found to be one order of magnitude 
greater than the coefficient c. 

Finally, when the correlation length is large, the above 
procedure leads to the approximate scaling form for the 
retarded function 


Xa/m(^5 ^ 


2 ^ 1 
Usiol + q^+^^ 


(39) 


the microscopic length scale 


^ ^_ 1 _ dx°f^{ci,iv) 

° Xa/m(q = 0>**^ = 0) 9q 

the mean field U for a phase transition 

Xa/m(q = 0 ,W = 0 )’ 


(41) 

(42) 


the deviation from the mean field U, 


6U = Umf-Us, 


(43) 


and 


J_ ^ _1_ 9x°afra 

ro ^0Xa/m(q = 0,^ = 0) duj 


(44) 


with xTfm imaginary part of the retarded suscepti¬ 
bility. 

For practical calculations, it is convenient to define the 
spin correlation lengths for the ferromagnetic and anti¬ 
ferromagnetic channels as 


^ X}„(q = 0 ,w = 0 ) 
xjm(^ = o,iu = oy 

^Xlfmi^ = 0,W = 0) 

xlf^ici = o,iu = oy 


(45) 

(46) 


where the correlation length is given by 


Indeed, using the scaling form Eq. (39), the above defi¬ 
nition corresponds to 





(40) 


Us^ox}mici = 0,w = 0) 

^ C Umf 

Co u, ■ 


(47) 


In these equations we have used the following definitions: 


( 48 ) 
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Since Umf/Ug ^ 1 when ^ is large, the two definitions of 
correlation lengths essentially agree in that limit. 

Although similar definitions of correlation lengths can 
be adopted in the charge channel, these lengths never 
become large so they are not really useful. 

We end with a note on the critical exponents. TPSC 
gives us a good estimate of the zero-temperature critical 
value of [/, although the exponents usually take values 
associated with the spherical model^^. Accurate values 
of exponents are usually found with the renormalization 
group approach. This is complementary to our approach 
since the latter methods do not give nonuniversal num¬ 
bers such as the critical U. For graphene, the universal¬ 
ity class is that of the Gross-Neveu model^^’^^ with 0.88 
as the value of the correlation length exponent to lead¬ 
ing order in e. Instead, we have the value 1, as follows 
from Eq. (40). From the scaling form (39), we see that 
the dynamical critical exponent defined by uj ~ ^ is 

z = 1. Lorentz invariance suggests that Fq equals the 
Fermi velocity vf^ while a better formula for the q and uj 
dependence in the denominator of Eq. (??) would prob¬ 
ably replace ^ by — {uj/vp^^ Further details 
appear in Appendix C. 


V. NUMERICAL PROCEDURE 

We first evaluate the noninteracting susceptibility 
(Lindhard function) x^(g) in Eq. (28). We then take a 
guess for to initialize the irreducible spin ver¬ 

tex Us^ Eq. (23). Using Eq. (26), we compute which, 
when substituted in the spin sum-rule, Eq. (29), allows 
us to update the variables {ua^Ua]) and since 

we know the filling (nc^cr) = 0.5 on the right-hand side. 
We repeat the procedure till we obtain self-consistent so¬ 
lutions for {ua^Uai) and and thereby obtain the 

irreducible spin vertex Us- 

A C++ code was written to calculate the noninteract¬ 
ing susceptibilities Xaa^ Xab^ Xba' FFTs are used in 
computations to exploit the convolutions in the defini¬ 
tions of the susceptibilities^^. First, the susceptibilities 
are computed in the position-imaginary time represen¬ 
tation where the convolution is just a product. EFT in 
the position space and a combination of cubic splines 
and EFT in the imaginary time space are implemented 
to obtain the final result in the momentum-bosonic Mat- 
subara frequency representation. The real (momentum) 
space grid is A x A, where N = 50, 100 and 200 were 
taken. Since the noninteracting susceptibility obeys, 

= = (49) 

Q 

we fixed the optimum value for the number of Matsubara 
frequencies by requiring that the above be satisfied to 
1% accuracy. Accordingly, the range of imaginary time 
from 0 to /d was divided into Np = 2 slices. 
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FIG. 2. (Color online) Plot of as a function of U for 

given temperatures(A = 100). The temperature dependence 
is very small, as can be checked from the inset. 


Further comments on finite-size effects and computa¬ 
tional procedure may be found at the end of Appendix A. 



FIG. 3. (Color online) Plot of the irreducible vertex for the 
spin Us as a function of U for given temperatures (N = 100). 
The temperature dependence is extremely small as can be 
seen from the inset. 


VI. RESULTS AND DISCUSSION 

Following the numerical procedure detailed above, 
we computed the double occupancy {ria^riai) self- 
consistently. This allowed us to obtain the TPSC spin 
susceptibility, Eq. (26), as well as the correlation length 
in the antiferromagnetic channel. 

Double Occupancy Due to bipartite symme¬ 
try, {ua^Uai) = which we define as Fig¬ 
ure 2 shows plotted as a function of interaction 

U for given temperatures. In the noninteracting case 

























[/ = 0, double occupancy factors into a product of the 
occupations for up and down electrons. At half-filling, 
(%,^) = 0-5, so that (%%) = {n^){ni) = 0.25 for U = 0. 
As U increases, the energy cost for two electrons occupy¬ 
ing a single site increases, thereby leading to a decreasing 
value of 

Spin vertex Us- Figure 3 shows the spin vertex Us as a 
function of the interaction U for given temperatures. 

For very small values of t/, Ug is almost equal to U. 
As U increases, Ug becomes less than U and it shows 
a tendency to saturate to a constant value. This is 
a result of Kanamori-Brueckner^^’^^’^^ screening: the 
physics reflects the fact that, as U increases, the two- 
body wave-function becomes smaller when electrons are 
on the same site to reduce the probability of double oc¬ 
cupancy, thereby decreasing the value of the effective on¬ 
site interaction. The maximum energy this can cost is 
the bandwidth so that at large values oiU^Ug saturates 
to a value of the order of the bandwidth. 

Correlation lengths for spin and eharge suseeptibili- 
ties. With the irreducible spin and charge vertices, we 
can calculate the spin and charge susceptibilities using 
the particle-hole Bethe-Salpeter equations (26) and (27). 
From the definitions Eqs. (35) and (36) we obtain the 
spin susceptibilities in the ferromagnetic and antiferro¬ 
magnetic channels and deduce the correlation lengths in 
the respective channels from Eqs. (45) and (46). 

The charge correlation lengths are not physically rele¬ 
vant since the irreducible charge vertex is generally larger 
than U and suppresses the charge susceptibility com¬ 
pared with its noninteracting value, meaning that the 
correlation lengths are always small and ill-defined. 



FIG. 4. (Color online) Semi-logarithmic plot of the spin cor¬ 
relation length in the antiferromagnetic channel ^ as a func¬ 
tion of U for various temperatures. (N = 100) 


Figure 4 shows the variation of spin correlation length 
in the antiferromagnetic channel ^ as a function of U 
for various temperatures. We first obtain the ratio of 
the interacting susceptibility to the noninteracting sus¬ 
ceptibility in the antiferromagnetic channel Ca/m ^^ing 



FIG. 5. (Color online) Plot of the ferromagnetic cor¬ 
relation length Cfm ^ function of U for given values of 
temperature. (A = 100). 


Eq. (46) and multiply it by the microscopic length 
(Eq. (48)) to obtain the correlation length f in units of 
the lattice spacing. The figure clearly indicates that the 
spin susceptibility in the antiferromagnetic channel in¬ 
creases steadily with increasing U and with decreasing 
T, as expected, with a clear tendency to diverge at suf¬ 
ficiently large U and low T. The quantitative accuracy 
of the results cannot be trusted for correlation lengths 
smaller than unity or larger than about half the system 
size. 

Eigure 5 shows the plots for the ferromagnetic corre¬ 
lation length estimated from the ratio of the inter¬ 
acting susceptibility to the noninteracting susceptibility, 
Eq. (45), as a function of U for various temperatures for 
N = 100. We can see that the ratio decreases as temper¬ 
ature decreases. Thus in the ferromagnetic channel, the 
correlation length never becomes larger than the lattice 
spacing and hence we do not focus on that case. 

^0 CLnd Fq We numerically determine and Fq for 
the spin susceptibility in the antiferromagnetic channel 
from the relations given in Eqs. (41) and (44) respec¬ 
tively. Eigure 6 shows the resulting temperature depen¬ 
dence of Fq and of ^o- The microscopic length scale 
is almost T independent and of order 0.25. Fq converges 
to the expected value vp = ^ only at very low T and 
for large values of A ^ 2000. Further discussion and 
numerical estimates appear in Appendix C. 

Crossover temperature and Uc- For U > Uc there is an- 
tiferromagnetism at T = 0. We expect then that, when 
U > Uc^ below a crossover temperature Tx, the anti¬ 
ferromagnetic correlation length becomes so large that 
one enters the renormalized-classical regime where the 
characteristic spin fluctuation frequency ujgf is less than 
temperature. And indeed, since the scaling form (39) 
implies that ujgf ^ and ^ increases faster than T 
at sufficiently low T when U is larger than Uc (Fig. 8), 
this implies that for U > Uc there is necessarily a tem¬ 
perature below which the condition ujgf < T is real- 
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FIG. 7. (Color online) Plots of the crossover temperature 
as a function of interaction strength for given values of N. 
The crossover temperature was determined from ^ = vf/T 
scanning U at hxed T, and scanning T at hxed U to obtain 
an estimate of the error in finding the intersection ^(U,T) = 
vfIT. 


FIG. 6. (Color online) Plot of (dashed blue line with 
circles) and of Fo (dot-dashed green line with squares) as a 
function of InT for the range T = 0.001 to T = 0.1. Each 
quantity has its own vertical axis: left for and right for Fq. 


ized. In the more standard case where the dynamical 
critical exponent satisfies 2 : = 2, a pseudogap in the 
single-particle density of states appears at a temperature 
smaller than that where Ugf ^ T. At that temperature, 
^ becomes larger than the thermal de Broglie wavelength 
vp/T (with vp the Fermi velocity).The question 
of the appearance of a pseudogap in the present case re¬ 
mains to be investigated, but it is expected as a precursor 
since there is a real gap in the antiferromagnetic state. 
The pseudogap should appear basically when we enter 
the renormalized classical regime since here frequency 
and wavevector scale in the same way. 

We thus define the crossover temperature to the renor¬ 
malized classical regime by the condition ^ = vp/T, with 
vp 8it the Dirac point. In order to extract the crossover 
temperature for a fixed value of [/, we plot the corre¬ 
lation length as a function of temperature and pick the 
value of temperature {Tx) where this plot intersects the 
plot of vpjT as a function of temperature. Similarly, 
for a fixed value of T, we can pick the value of inter¬ 
action (Ux) where the correlation length exceeds vp/T. 
Figure 7 shows the plots of crossover temperature as a 
function of interaction determined using both approaches 



0 

II 

0 

0 

II 

N = 200 

T vs. Ux 

linear 

3.8 

3.794 

3.793 

quadratic 

3.825 

3.806 

3.808 

Tx vs. U 

linear 

3.779 

3.775 

3.775 

quadratic 

3.809 

3.795 

3.789 


TABLE I. Values of the critical interaction strength Uc ob¬ 
tained from the linear and quadratic extrapolation of the 
crossover plots in Eig. 7 for various values of N. 


detailed above, for N = 50, 100, and 200. By quadratic 
and linear extrapolations of the curves to zero tempera¬ 
ture, one obtains the results for Uc that appear in Table 

I. 

Critical exponent z and an alternate determination of 
Uc- We can find the critical value Uc using another ap¬ 
proach. This approach lets us estimate the dynamical 
critical exponent 2 : also. In Fig. 8 we plot In ^ as a func¬ 
tion of InT for N = 100 and 200 where the correlation 
length is sufficiently small that finite-size errors are not 
important (except far from Uc). For U < Uc^ In^ satu¬ 
rates at low temperatures, while for U > Uc, Inf diverges 
and finally, at Uc, f has a pure power law behavior. In 
order to determine Uc, we fit In^ versus InT for vari¬ 
ous values of U with straight lines. The value of U that 
gives the best fit is taken as Uc. It is the slope of In^ 
vs InT that gives us the numerical estimate of the dy¬ 
namical critical exponent Despite the fact that we are 
not in the asymptotic regime for ^0 and Fq, the value 
so obtained is 2 : = 1.00 for Uc = 3.8 ± 0.005. For high 
temperatures, all curves have the same slope as the case 
U = Uc. For A" = 50, where we saw finite-size effects 
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in Fig. 7, the largest correlation length is close to N/2 
at the smallest temperature for U = Uc, invalidating the 
estimate. Indeed, in that case Uc = 3.85 ± 0.005 but 
2 : = 0.87, which is clearly incorrect. 

Taking the average value of Uc obtained for = 100 
and N = 200 in Table I and estimating the error from the 
range of values obtained, we find that Uc = 3.79 ± 0.01, 
consistent with the result obtained from the estimate of 
the previous paragraph with z = 1. 


renormalized classical regime and where a pseudogap is 
expected to open up. 

Generalized extensions of TPSC to multiband cases of 
the type presented here and in Ref. 58 have the poten¬ 
tial to open the study of interacting systems, and to im¬ 
prove realistic materials calculations. In the latter case, 
TPSC offers the possibility to include long wave length 
spin fluctuations in addition to long wave length charge 
fluctuations already present in these approaches. 



FIG. 8. (Color online) Plots of In^ as a function of InT for 
various values of U {N = 100 and 200). The straight dashed 
magenta line corresponds to a pure power law, z = 1, hence 
to the value Uc for the quantum critical point. 


VII. CONCLUSION 

The nonperturbative TPSC theory has been extended 
to a multi-band case, namely the Hubbard model on the 
honeycomb lattice. In TPSC, valid from weak to in¬ 
termediate coupling, charge and spin irreducible inter¬ 
actions are determined self-consistently in such a way 
that conservation laws and the Pauli principle are satis¬ 
fied. The Mermin-Wagner theorem is also automatically 
satisfied and the physics of Kanamori-Brueckner screen¬ 
ing that renormalizes the spin and charge irreducible 
vertices is taken into account. On the honeycomb lat¬ 
tice, nearest-neighbor antiferromagnetic fluctuations are 
dominant. The TPSC value of Uc/t for the quantum- 
critical semimetallic to antiferromagnetic transition is 
Uc/t = 3.79 ±0.01 consistent with^^ Uc/t = 3.869 ±0.013 
and^^ Uc/t = 3.78 obtained from the large scale quantum 
Monte Carlo calculations and also consistent with the 
functional renormalization group^^’^^ Uc/t = 3.8. These 
results rule out the existence of a spin liquid phase in the 
ground state of the graphene Hubbard model at interme¬ 
diate couplings since estimates for the Mott transition 
yield a Umou larger than Uc- We have also estimated 
the crossover line in the T-U plane where one enters the 


ACKNOWLEDGMENTS 

We are grateful to Dominic Bergeron and Wei Wu for 
illuminating discussions. S.A. would like to thank P. 
Mangalapandi for his expert advice on parallel computa¬ 
tion. This work was supported by the Natural Sciences 
and Engineering Research Council of Canada (NSERC), 
and by the Tier I Canada Research Chair Program (A.- 
M.S.T.). 


Appendix A: NONINTERACTING 
SUSCEPTIBILITIES 


Consider the spin susceptibility 


Xap{h2) = {T,Sl{l)Sl{2)) (Al) 

where = a^b and S'^(l) = ~ %,«(!)• Eval¬ 

uating this expression in terms of noninteracting Green 
functions in Eq. (28), we find 


2 _^ 2 

= -Jp E 4 [^a/3(k,q,w)] 

k,a,/3 


where is a Matsubara frequency and 

^a/3(k,q,2Z^) = --- 

and (a, /3 = ±, with the Eermi function 


1 


niE^) ^E^/T ^ ^ 


and eigenenergies 

E^ = a\f{k)\ 


(A2) 


(A3) 


(A4) 


(A5) 


= (ay /3 ± 2 cos ki 2 cos ± 2 cos(—— /C 2 ). 


(A6) 


/(k) is defined in Eq. (4), T is the temperature, while 
ki and k 2 are the components of the momentum vector 
on the two unit lattice vectors ai and a. 2 . The other 
components of the susceptibility tensor are given by 
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Xa6(q>*'^) = - [al3Ma0{Kfl,ii')] 

k,a,/3 


where = |^|^|| and 


xL(q> [Q;/3M„,^(k, q, w)]. 

k 


(A7) 


(A8) 


The relation between Xba Xab (33) is thus satisfied. 


The equivalence of the two sublattices implies Xaa — Xbb 
but in general not Xab ~ Xba because of the chiral nature 
of the Dirac points. 

Instead of using FFT’s, one can first perform the Mat- 
subara frequency sum exactly and then sum over wave 
vectors. In that case, there are certain points in the Bril- 
louin zone where the Mq,^ have the form 0 / 0 . At these 
points one must take the limits and use rHospital’s rule. 
For example at q = 0 


q^O 


^ , n{E+)-n{E+J 

0) = lim X \ 

<,^0 E+-E+^, 

(A9) 

q^O dE^q 

(AlO) 

= -/3n{E+)[l-n{E+)]. 

(All) 


Similarly = 0/0 at k = (27r/3,27r/3) and q = 

(27r/3, 27r/3) so the same solution applies. However, this 
procedure means that the Dirac points introduce large 
finite-size effects in the temperature dependence. Choos¬ 
ing a grid that is regular but avoids the Dirac points (for 


I 

example N = 100 instead of A" = 90) minimizes finite- 
size effects. 


Appendix B: ANTIFERROMAGNETIC 
SUSCEPTIBILITY 

In the presence of interactions the spin susceptibility 
is given by the matrix equation 


A 


V = (I- 


-1 


(Bl) 


Defining the antiferromagnetic susceptibility by 

Xa/m(q, w) = Xaa(q, w) + \/xo6(q, w)x6o(q, w) (B2) 

allows US to find a simple scalar equation that reduces 
to Xafm = Xaa + Xab whcu Xab IS real. The combination 
Xaa + Xab docs uot Satisfy a simple scalar equation in the 
general case. 

The algebra that follows proves these assertions. Ex¬ 
panding the matrix equation we find 


= 


J_( 

^ 1 - Xxtb 

%X°ft ] 

( Xaa 

Xab 

det ' 

\ 2 ^ba 

1 - XxL ) 

V X°a 

Xbb 



( i^-XxWaa+XxlbXl 

Xab 

det ' 

( Xba 

fxLx°b + (i-^xMb, 


where det, that stands for the determinant, can be expanded as 

det = (1 - YX6b)(l - yXaa) - ^xlbxta 


t 2 ^Xbb ”1” Xaa) ^ iXabXba XaaXbb) 

= - y(\/xLX° 6 + /XabXL)) (l - “ /XabXba 


With Xaa — Xbb simplify the determinant 


det = ( 1 - Y(Xaa 


+ 


\JXabXba)) (^1 - Y^Xaa “ ]/XlbXba) 


(B3) 

(B4) 

(B5) 

(B6) 

(B7) 

(B8) 

(B9) 
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and the antiferromagnetic susceptibility 


^afm Xaa 's/Xab^bi 


C 


(1 - ^Xbb)Xaa + ^xlbXba + VXgbXbg 

1 - ^{XL + VXabXba)) (l - ^(XL - 'Jx'ibXbc 

.0 ^ ^ 0 _ U^((^0 _ ^0 0 

V ^ah^ha 2 VVAaa/ ^ab^ba) 

^ /n^Tn“^^ (a /tj 


_ ' '^ab'^ba _ 2 vv-^aa/ _ '^ab'^ba/ _ 

(i - ^{xL + VxlbXba)) (i - ^(xL - VxlbXbc 

(Xaa + VXabXba) “ ^iXaa “ VXabXD) 

1 - ^ixL + VXabXba)) (l - ^ixL - VXabXbc 

..0 , /_0 _0 




yO 

Aaa 


VXabXb^ 

{^-^{xL + VXabXbai) 


that 


Xa/m(q >*0 = 


Xa/m(q^^O 
- ^Xa/™(q ,*0 


(B15) 


(BIO) 

(Bll) 

(B12) 

(B13) 

(B14) 


Appendix C: DERIVATIVES OF 
NONINTERACTING SUSCEPTIBILITIES IN 
THE DIRAC APPROXIMATION AND 
ESTIMATES FOR fo,ro 

1* Xafm Sind its derivatives 


At q = 0, the phase in Xab^ disappears and we are left with Xabi^ ~ O^iu) = Xbai^ ~ O^iu). Given that 

we are taking the positive square root, we also have ^/xab{0^ '^^)Xba{0^ = —Xab{0^ so that the retarded function 

is 


Xa/m(0,W + i(5) 


Xaa(0, UJ + i5)- x°b(0, w + iS) 


[■^+- q’ + ^-+ q’ + *'^)] 

k 

_Lw . ~^^(-5'k) + i 

■— . 2 -E'k + ^ + i5_ 


(Cl) 

(C 2 ) 

(C3) 


Only interband transitions contribute to Xo/m(0o^)- 

In the Dirac approximation, we evaluate separately the real and imaginary parts. Beginning with the latter, we 
find 


^^xlfm {0,oj + iS) = -TT^y^tanh {/3E+/2) (<5 {2E+ +w) -(l( 2 £;+ -w)). 


(C4) 


Transforming the sum into an integral, going to cylindrical coordinates, we have 


1 P dPk kdk 1 

J {2TTf ^ Jo 27r ^ v% Jo 


ede 

2it 


(C5) 
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where A^; is the energy cutoff. Assuming cj > 0, only the last S function contributes. Taking into account a factor of 
2 for the two Dirac points, we are left with 


2 eds 

ImXa/m(0,w + i(5) =7r^ / — tanh (/3£/2) <5 (2e - w). 

Jo 

In the zero temperature limit, 

2 r^E 

ImXa/m( 0 ,w + i(5) = ^ / 5{2e-w)ede 

Jo 

1 UJ 

^ ^ 2 ■ 

For the real part, we begin with 

Re xlfm (0> II tanh {^E+ /2) -T • 

Expanding around the two Dirac points as above gives 

^ 2 Ed£ Ae 

ReXa/m(0,w + i(5) =P;^ / — tanh(/?e/2) —— 2 - 

Working in the zero temperature limit, we find 


(C 6 ) 

(C7) 

(C 8 ) 

(C9) 

(CIO) 


Rex°/„(0,w + i(5) f 

Jo 


de (2ef 
V {2sf -uj^ 


I de 2 1 -T, f 

— / — + F—r / 
Vp Jo TT ^'F Jo 

I 


dE 1 


Assuming cj > 0 , we find 


^Af 


1 


F JO (2^) — 

2 ['2KeI^ 2 


2ujv‘^f 


dx- 


x‘^ — 1 


+ 2 ^ 2 , tt 

1 A 

-A, 


tanh 


-1 


2 A F 


FVp Av^fAe 


(Cll) 

(C 12 ) 

(C13) 

(C14) 

(C15) 


2. Estimates for and To 


We begin with the definition Eq. (44) of Eq and substitute the results just found Eqs. (C 8 ) and (C15) to find 


1 

f~0 


1 1 d 

^o^ex°fm (0,0) aw 

1 

I Av% _ F 


Im xlfm (0,w + i(5) 


Co ;^Af 4^oAf 


(C16) 

(C17) 


SO as expected Eq has units of velocity since with /i = 1, cutoff, then 

Af is {time)~^ . Taking Af = vfA with A = n/a the . 

44oAf . 4o 
Eo =-= 4 —ff 

F a 


(C18) 
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From the Lorentz invariance, we expect 

To = (C19) 


which, with a = 1, suggests that ^ 0.25. The result 
found numerically for in Fig. 6 is just slightly larger 
because band curvature means that Ke is a bit smaller 
than the estimate Ke = ffA. Similarly, Fq at low tem¬ 
peratures is numerically close to ve = \/3/2 in our units. 
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